Global attractors and extinction dynamics of cyclically competing species 
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Transitions to absorbing states are of fundamental importance in non-equilibrium physics as well 
as ecology. In ecology, absorbing states correspond to the extinction of species. We here study the 
spatial population dynamics of three cyclically interacting species. The interaction scheme comprises 
both direct competition between species as in the cyclic Lotka-Volterra model, and separated selec- 
tion and reproduction processes as in the May-Leonard model. We show that the dynamic processes 
leading to the transient maintenance of biodiversity are closely linked to attractors of the nonlin- 
ear dynamics for the overall species' concentrations. The characteristics of these global attractors 
change qualitatively at certain threshold values of the mobility, and depend on the relative strength 
of the different types of competition between species. They give information about the scaling of 
extinction times with the system size and thereby the stability of biodiversity. We define an effec- 
tive free energy as the negative logarithm of the probability to find the system in a specific global 
state before reaching one of the absorbing states. The global attractors then correspond to minima 
of this effective energy landscape and determine the most probable values for the species' global 
concentrations. As in equilibrium thermodynamics, qualitative changes in the effective free energy 
landscape indicate and characterize the underlying non-equilibrium phase transitions. We provide 
the complete phase diagrams for the population dynamics, and give a comprehensive analysis of the 
spatio-temporal dynamics and routes to extinction in the respective phases. 



Absorbing states play an important role in ecology, 
where they correspond to the extinction of species pp. 
While any stochastic system will eventually end up in 
one of its absorbing states, in nature, one finds a surpris- 
ingly long-term maintenance of biodiversity in ecosys- 
tems containing a broad variety of coexisting species. A 
structured environment in combination with cyclic com- 
petition between species was proposed to be a main facil- 
itator of biodiversity [2j |3] . Classical ecological examples 
for cyclic interactions comprise coral reef invertebrates 
[4], rodents in the high arctic tundra in Greenland [5], 
and cyclic competition between different mating strate- 
gies of lizards [6 j. However, long reproduction times and 
large spatial scales involved make it difficult to quanti- 
tatively analyze these ecological systems. To circumvent 
these problems, recent experimental studies have turned 
to microbial model systems comprising three genetically 
distinct strains of Escherichia coli which cyclically dom- 
inate each other like in the children's game rock-paper- 
scissors pJIE]. 

These experimental studies of microbial model systems 
have motivated a large body of theoretical literature ex- 
ploring the role of cyclic interactions in ecological sys- 
tems [TJ [3l [9H42] . Most of this work has focused on 
two paradigmatic examples of three-species models with 
cyclic interactions. In a first class of models, direct com- 
petition between two individuals leads to the immedi- 
ate replacement of the weaker species by the stronger 
one [3l [TTH181 [20H32]. This type of competition, where 
selection and reproduction are combined into a single 
process, is similar as in the classical two-species Lotka- 
Volterra model [43-45 j. The interaction scheme of this 
cyclic Lotka-Volterra model may be summarized by a set 
of chemical reactions between the three species A, £?, and 



C: 

AB AA, BC -+ BB, CA^CC. (1) 

In the second class of models, originally proposed by May 
and Leonard [19] , selection and reproduction are two sep- 
arate processes. An interaction between two individuals 
with different traits (strategies) leads to the death of the 
weaker species and thereby to empty spaces. Reproduc- 
tion then follows as a second process which recolonizes 
this empty space. The ensuing reaction scheme reads: 

AB^AQ, BC^B®, CA^rC%, (2a) 
A0 ->> AA , B$->BB, C$->CC. (2b) 

Both of these models exhibit absorbing states where 
all but one species have died out. Due to the inevitable 
demographic fluctuations in systems with a finite number 
of individuals these absorbing states will with certainty 
be reached if one just waits long enough. How long one 
has to wait strongly depends on the type of model and 
the ecological scenario under consideration. 

In well- mixed systems, the typical extinction time T 
was found to scale linearly with the population size N 
for the cyclic Lotka-Volterra model QU [H H2 US HZ| 
and logarithmically for the May-Leonard model [34 . The 
reason for the difference is the different nature of the 
phase space orbits characterizing the nonlinear dynam- 
ics of these two models pp. While the phase portrait of 
the cyclic Lotka-Volterra model exhibits neutrally sta- 
ble orbits, the May-Leonard model is characterized by 
heteroclinic orbits emerging from orbits which spiral out 
from an unstable reactive fixed point. For neutrally sta- 
ble orbits, the stochastic dynamics performs an unbiased 
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random walk which implies that T ex N. In contrast, un- 
stable orbits generate a drift of the trajectories in phase 
space towards the boundary such that the extinction pro- 
cess towards the absorbing states is exponentially accel- 
erated with T oc In N [TJ [T6j HE] . 

In spatially extended systems, the scaling of T with 
population size strongly depends on the degree of mix- 
ing. In particular, it has been shown for both models 
that there exists a mobility threshold below which ex- 
tinction times scale exponentially in the system size. For 
the May-Leonard model this has been attributed to the 
existence of spiral waves, which emerge as a result of the 
local nature of reactions and internal noise [33] [34[ [36] . 
Above a certain mobility the characteristic wave length of 
the spirals exceeds the system size, effectively rendering 
the dynamics well-mixed. In this regime, extinctions oc- 
curs rapidly. In the cyclic Lotka-Volterra model, spatial 
patterns are unstable as a result of an Eckhaus insta- 
bility [26]. However, below a mobility threshold biodi- 
versity is still maintained by strong spatial correlations. 
Further work has extended these findings to asymmetric 
reaction rates [14] [29] and more complex interaction net- 
works [32] [481- 50 . In a niche model it has been shown 
that interaction networks with a high connectivity and 
a hierarchical or cyclic interaction structure lead to in- 
creased diversity [51] [52] . For the May-Leonard and the 
cyclic Lotka-Volterra model it was found that spatially 
inhomogeneous reaction rates have only minor effects 
on the dynamics [27] [39] . For the classical two-species 
Lotka-Volterra model, analytical studies have been per- 
formed to understand the underlying mechanism lead- 
ing to the stabilizing correlations [53] [54]. These studies 
argue that the stabilization can by understood by the 
desynchronization of diffusively coupled oscillators. The 
desynchronization is a result of the combined effect of 
noise, migration and the dependence of the oscillations' 
frequency upon their amplitude. 

For the one-dimensional May-Leonard model the dy- 
namics leading to extinction has been studied in greater 
detail. If the individuals diffuse only little or do not dif- 
fuse at all, coarse-graining of species' domains has been 
identified as the dominant dynamical process leading to 
extinction [20[ |23[ [24] [27] . With increasing diffusion con- 
stant other types of collective excitations become impor- 
tant [38 . The dynamics to extinction is then surpris- 
ingly rich, comprising rapid extinction, global oscillatory 
behavior, and traveling waves. The latter involve oscillat- 
ing overall densities, i.e. the domain sizes for the differ- 
ent species change periodically. The statistical weights of 
these dynamical regimes change qualitatively at thresh- 
old values of the mobility and the system size. Taken 
together, it has turned out that the dynamics in the one- 
dimensional May-Leonard model is highly complex, much 
more than one would naively anticipate. 

In this manuscript we extend these studies to two- 
dimensional models with cyclic competition between 
species. Specifically, we study a generic model comprising 
both direct competition between species as in the cyclic 



Lotka-Volterra model, and separated selection and re- 
production processes as in the May-Leonard model. Our 
goal is to identify and characterize the dynamic processes 
which are responsible for the transient maintenance of 
biodiversity and which finally lead to the extinction of 
all but one species. In particular, we are interested in 
how factors like species mobility and the relative strength 
of the different competition types govern the complex 
spatio-temporal dynamics of the system. Employing ex- 
tensive numerical simulations, we show that the dom- 
inant dynamic processes responsible for the transient 
maintenance of biodiversity correspond to attractors of 
the global dynamics, i.e. the overall density of species 
in the system. The characteristic features of these at- 
tractors give information about the scaling of extinction 
times with the system size and thereby the stability of 
biodiversity. Importantly, the attractors change quali- 
tatively at certain threshold values of the mobility and 
the relative strength of the different competition types. 
The phase transitions at these threshold values corre- 
spond to abrupt changes of the scaling of the extinction 
time T with population size N. These global attractors 
can be envisioned as minima in an effective free energy 
landscape. As their counterparts from equilibrium ther- 
modynamics, they give valuable information about the 
physics underlying the observed transitions and thereby 
give insight into the mechanisms leading to the stability 
of ecosystems. Our numerical studies are complemented 
by scaling arguments based on properties of the complex 
Ginzburg-Landau equation [26, 33, 34, 36, [55]. 

I. A GENERIC MODEL OF CYCLICALLY 
INTERACTING SPECIES 

We consider a spatially extended population consisting 
of three distinct species A, B and C that compete with 
each other cyclically in two different ways: either by im- 
mediately replacing the competitor by an individual of 
its own kind, or by killing the inferior species and cre- 
ating an empty site 0. In addition, individuals may also 
reproduce if empty spaces are available. These processes 
are summarized by the following reaction scheme: 

AB A A0, AH) AAA, AB^AA, (3a) 
BC^B®, Bty-^BB, BC^BB, (3b) 
CA^>C0, C0 4CC, CA^CC. (3c) 

The reaction rules ([3| describe two competing types of 
selection processes: On the one hand, with rates a and 
/i, selection and reproduction are separate processes. Se- 
lection produces empty sites which are in turn required 
for reproduction. An empty space is not necessarily occu- 
pied by the individual who produced it. We refer to these 
processes as May-Leonard processes. On the other hand, 
Lotka-Volterra processes, with a rate couple selection 
and reproduction: success in competition directly trans- 
lates into reproduction. In the following, when we use 





the term Lotka-Volterra process, this will always imply 
that the reactions are cyclic. There are two limiting cases 
which correspond to well established models: for v — » 
and for \i = a = we recover the May-Leonard model 
[T9] and a three species model with cyclic interactions of 
Lotka-Volterra type [9] [43] [44] , respectively. 



A. Stochastic lattice gas model 

We consider a two-dimensional square lattice and em- 
ploy periodic boundary conditions [56]. The linear di- 
mension of the lattice is taken as the basic length unit 
such that the lattice constant a = 1/L with L the num- 
ber of lattice sites along each axis. At each site a fixed 
number M of individuals (A, 5, C or empty spaces 0) 
are located. M may be viewed as the carrying capacity 
of a lattice site. In addition, individuals are also able to 



move on the lattice. While the reactions, Eqs. (3a)- (3c 
are assumed to occur on the same lattice site, the indi- 
viduals' mobility is implemented as an exchange process 
at a rate e between neighboring sites, XY — >• FX, where 
X and Y denote species A, B and C or empty spaces 
0. Macroscopically the nearest neighbor exchange pro- 
cess leads to diffusion with an effective diffusion constant 
D = eL~ 2 /2 [33 El [36]. As two particles are involved 
in migration, it also induces some additional nonlinear 
reaction terms, which we neglect here [57|l58]. 

We performed extensive simulations of the ensuing 
stochastic particle dynamics employing a sequential up- 
dating algorithm: At each simulation step an individual 
is chosen at random. It then either reacts with another 
also randomly chosen individual from the same site, or is 
exchanged with an individual of a neighboring stack; each 
stochastic event occurs with probabilities corresponding 
to the respective reaction rates /i, a, v and e. Typi- 
cal snapshots of the stochastic simulations for the May- 
Leonard model are shown in Fig. [T] 

The effective size of the system is N = M ■ L 2 . If 
M and L are large enough, the strength of fluctuations is 
proportional to 1/ y/N [59 . The simulation results there- 
fore do not depend on the specific choice of M and L, as 
long as both are not too small and the net system size is 
kept constant. In particular, the lattice spacing a = L~ x 
must be much smaller than the correlation length £. 

Different reaction rates for the species should not limit 
the validity of our results, as long as the differences be- 
tween the species are not too large. It has recently been 
shown that small asymmetries in the reaction rates do not 
alter the dynamics [35]. A general discussion of species- 
dependent reactions rates is given in Refs. [27] EH]. In 
the following we will also set n = a. While the relation 
between the selection and reproduction rates in the May- 
Leonard model affects certain properties of the dynamics 
(such as the wavelength and velocity of spiral waves), 
qualitatively the results remain the same [36] . This view 
is supported by some sample runs that we carried out 
for different values of fi/cr. It is, however, important to 




FIG. 1. (Color online) Typical spatial patterns in the May- 
Leonard model. Color (gray scale) denotes the concentra- 
tion of species A, B and C, with red (medium gray) signify- 
ing a site dominated by species A, green (light gray) a site 
dominated by species B, and blue (dark gray) being a site 
dominated by species C. (a) For large diffusion constants 
(D = 1.5 • 10 -3 ), the dynamics shows global oscillations with 
periodic switching between states with mainly one species 
present, (b) For intermediate values of the diffusion constant 
(D = 5-10 -4 ), we observe planar traveling waves. Here, two of 
the domains take a characteristic domain size dictated by the 
diffusion constant. The third domain then occupies the rest 
of the system, (c), (d) For even smaller diffusion constants 
(D = 3 • 10 -4 , and D = 6 • 10 -5 ), pairs of rotating spirals 
appear. The vertices of these spirals move very slowly on a 
time scale much larger than the time scale corresponding to 
the propagation speed of their arms. The wave length of the 
spirals decreases when the diffusion constant is reduced. The 
system size for all snapshots was L = 60, and the carrying 
capacity for each site was chosen to be M = 8. 



note that our results are not valid for extreme choices 
of the rate, corresponding, for example, to two species 
predator prey models [45] [60j [61] . In all simulations the 
initial condition was chosen to be a randomly populated 
lattice with average concentrations corresponding to the 
reactive fixed point of the well mixed model. 

We fixed the time scale by setting /i = a = 1 for the 
May-Leonard and v = 1 for the Lotka-Volterra limit. The 
diffusion constant D then gives the mean square displace- 
ment of an average particle between two reactions. As an 
example, with the system size as the unit length a value 
of D = 10 -3 implies that a particle covers an area of one 
thousandth of the system size between two succeeding re- 
actions. We study a regime, where the correlation length 
is much larger then the lattice spacing. As a result, the 
lattice spacing is irrelevant as a length scale. The only 
relevant quantity is the ratio of the characteristic length 
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of spatial patterns given by D and the system size. 



B. Well-mixed limit and invariant manifolds 

For large populations, intrinsic fluctuations are negli- 
gible in our stochastic lattice gas model. If in addition 
every individual can interact with every other individual 
with equal probability, i. e. for a well- mixed system, the 
dynamics is aptly described by deterministic rate equa- 
tions for the densities a, b and c of the species A, B and 
C, respectively: 



d t a = 
d t c 



-aba ■ 
-acb - 



- [10,(1 

-pb(l - 
■ pc(l - 



-p) 
■ P )- 

p)- 



- va(b - 
■ vb(c - 
vc(a - 



-c), 

a) , 

b) . 



(4a) 
(4b) 
(4c) 



Here, p = a + b + c is the total density of species and 
< a, 6, c, p < 1. Equation Q exhibits three absorb- 
ing fixed points: (1,0,0), (0,1,0), and (0,0,1). They 
correspond to the extinction of two of the three species. 
Another fixed point at (0,0,0) corresponds to the ex- 
tinction of all three species. However, this fixed point 
cannot be reached by the stochastic dynamics from the 
initial conditions we study here. In addition, there is a 
reactive fixed point 



(a*,b*,c*) = 



P 



(1,1,1), 



(5) 



at which all three species coexist. The dynamics in the 
vicinity of the reactive fixed point can be studied by lin- 
earizing Eq. Q around (a*,6*,c*) and by determining 
the eigenvalues of the corresponding Jacobian. We find 
that the dynamics close to the reactive fixed point is char- 
acterized by an attractive eigendirection with a negative 
eigenvalue k,q = —p and two further eigendirections with 
eigenvalues 



K± 



1 p 



2 3p 



(6) 



Therefore, the eigenvectors corresponding to k± span, 
to linear order, an invariant manifold onto which the dy- 
namics relaxes exponentially fast. To obtain an approx- 
imation for the invariant manifold, valid to second or- 
der in the concentrations, we follow the steps given in 
Ref. [36 . We first transform to a new reference frame 
whose origin is the unstable fixed point, (xa,%b,%c) — 
(a — a*, b — 6*, c — c*). Further, we choose the eigendi- 
rections of the fixed point as basis vectors for our new 
reference frame. To this end, we employ a rotation of the 
coordinate system, 



/ o -Vs ' 

-12 -1 

V i i i 



(7) 



The stable eigendirection corresponding to is then 
given by the yc~ direct ion, while ua and ys span the in- 
variant manifold to linear order. We parametrize the 



invariant manifold by yc = G(yA-,yB)- Using the ansatz 
G(yA, Vb) ~ y\ J ry 2 B an d determining the proportionality 
constant such that 



dtG(y A (t),y B (t)) = dy A G-d t yA+dy B G-d t yB = d t yc 

we find that G^a^b) is, to second order, given by 
a 3p - 



G(y A ,VB) = 



4/i 3p + 2a 



(va + Vb) - 



yc=G 



(8) 



This equation is valid only for p ^ 0. In the limit 
of a cyclic Lotka-Volterra model, p = = a, we find 
G fl= o j(J= o(yA^yB) — and the invariant manifold is 
given by the unit simplex defined by a + 6 + c = 1. 
This result can also be directly inferred from the Lotka- 
Volterra reactions, which preserve the total density and 
thereby lead to dynamics on an invariant manifold given 
by a(t) + b(t) + c(t) = 1. 

The rate equations in the new reference frame read: 



dtVA = —(2v + cr) (y\ 



Vb) 



(9a) 



V3 
2 

VA{pcr 



(2v + o)y B 
(3/i4 



3p + a 

o) Wvb - 



-yc 

(6p + a)y c ]} 



2(3/i + a) 
y/3(2v + a)y B [p + (3p + a)y c ] 



2(3/i + a) 



dtVB 



-r 2 v A 



vl) 



(9b) 



V3(2 



v-\-cr 
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Vb 



2(3p + a) 



[pa 



p 



Vb 



3p + ' 
(3p + a)(6p 



d t yc = -mc - (3/i + o)y\ 



(y 2 A 



rye 

cr)y c 

yl). 



(9c) 



What is the simplest differential equation that captures 
the essential features of the rate equations Q? Such 
a differential equation is called normal form, and is ob- 
tained by a nonlinear transformation which eliminates 
the quadratic terms. Following the steps in Ref. [36 , one 
makes a quadratic ansatz for the transformation and de- 
termines the coefficients canceling the quadratic terms. 
We find that the transformation is given by 



ZA 



zb 



va 
vb 

with prefactors 

OL\ - 
OL2 - 



(V3y 



(OL 2 2 

«i [y Va 
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(10a) 
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5 



Upon introducing a complex variable z = za + i^s and 
neglecting terms of order 0(z 4 ), the dynamics can finally 
be written in the form 



rules: 



where 



c 2 



d t z = (ci - io;)z - c 2 (l - ic 3 ^ 



\/3/i(2^ + cr) 

1 /i<7 

2 3/i + a ' 
cr(3/i + cr) 



\z\ 2 z. 



C3 



56/i(3/i + 2a) 

(J 2 (48/i + llcr) + 3^(60/i + 13cr)0 + a) 
a 2 + f v{v + cr) " 

1 (3/i + cr)v / 3(2z^ + cr) 
C2 56/i(3/i + 2a) 
cr 2 (18/i + 5a) + 9z/(6/i + a)0 + a) 

r2 I 27 



(12) 

(13a) 
(13b) 
(13c) 

(13d) 



While the limiting case of a May-Leonard model is found 
by simply setting v = 0, the cyclic Lotka-Volterra model 
is recovered by first taking the limit a — >• and then 
li 0. Other ways for performing this limit are pos- 
sible. However, taking the limit in a first, ensures that 
we obtain the established cyclic Lotka-Volterra model, 
which does not comprise empty sites: a + b + c = 1 [62] . 



C. Spatially extended continuum model 

In a continuum formulation, the nearest neighbor ex- 
change process macroscopically leads to diffusion with a 
diffusion constant D = eL~ 2 /2. The ensuing diffusion- 
reaction equations are simply obtained from the rate 
equations Q by supplementing them with diffusion 
terms DV 2 a, DV 2 b, and L>V 2 c, respectively j26 ] l33 | l34 l 
[36] . Then, upon applying the above transformations to 
the diffusion terms one obtains 

dtz = D V 2 z — i (Vz*) 2 + reaction terms, (14) 

where we have neglected gradient terms of order 
(D((\7z) 3 ). We expect that the dynamics is dominated 
by the long wavelength modes, and therefore only keep 
the leading order gradient term, leading to normal diffu- 
sion in the complex concentration z. This finally leads 
to the complex Ginzburg-Landau equation, 

d t z = DV 2 z + (ci - \uj)z - c 2 (l - ic 3 )\z\ 2 z , (15) 

a paradigmatic equation in nonlinear dynamics [55] . 



II. THE MAY-LEONARD LIMIT 

The May-Leonard model, obtained in the limit v — >• 0, 
is characterized by the following reduced set of reaction 



AB A A0 , 
A0 A AA , 



BC^B®, CA^C®, (16a) 
B$ A BB , C0 4CC. (16b) 



For large systems in the well- mixed limit, the dynamics 
is described by the May-Leonard equations [19], 



d t a — —aac + fia(l — p) , 
d t b = —aba + /x6(l — p) , 
dtc = —crcb + /ic(l — p) . 



(17a) 
(17b) 
(17c) 



The nonlinear dynamics of these equations is charac- 
terized by the same types of fixed points and invari- 
ant manifold as the general model Q. The reactive 
fixed point (a*,6*,c*) is globally unstable, as manifested 
by the existence of the Lyapunov function C = abc/p 3 . 
When starting in the vicinity of the unstable fixed point, 
the trajectories spiral outward on the invariant mani- 
fold and form heteroclinic cycles, approaching the bound- 
ary of the phase space and the absorbing states without 
ever reaching them [19]. However, intrinsic noise due 
to the stochastic nature of the interactions, and spatial 
structure drastically alter the observed behavior. While 
in well-mixed systems stochastic fluctuations drive the 
system into one of the absorbing states within a short 
time proportional to the logarithm of the system size 
HD CE2 El US CE3 CEH1 E2], spatial structures may effec- 
tively delay extinction by orders of magnitude [3] [34] . 

Similar to the previously studied one-dimensional 
case [38 , the two-dimensional, stochastic May-Leonard 
model exhibits distinct dynamical regimes as a function 
of the diffusion constant D (Fig.[T]). From our simulations 
we find the following phenomenology: For large values of 
D, we observe that the system (after some initial tran- 
sient) is first almost entirely taken over by one species, 
but with a few individuals of a second species surviv- 
ing, which dominates over the more abundant species. 
This second species will then slowly take over the sys- 
tem and thereby lead to a dynamical behavior that is 
reminiscent of the heteroclinic orbits of the determinis- 
tic, well-mixed system, where the global dynamics ap- 
proaches the boundary of the invariant manifold. In this 
regime, spatial patterns are of minor importance and the 
dynamics can be understood in terms of a quadratic coag- 
ulation process as outlined in Ref. [38 . With decreasing 
diffusion constant we observe planar waves of cyclically 
aligned uniform domains as well as rotating spiral waves 
[Fig. [TJb)-(d)]. In planar waves the overall concentra- 
tions may be constant, corresponding to stable domain 
borders, or change periodically, as a result of "tunneling" 
events in the leading edges of the fronts. The leading 
edges of the fronts may reach into second next domains, 
i.e. there is a finite probability for particles to pene- 
trate domains of prey via "tunneling" events [38]. As a 
consequence domain sizes oscillate periodically between 
characteristic length scales, thereby leading to oscillat- 
ing overall densities. The dynamics of rotating spirals 



D=9.540 
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FIG. 2. (Color online) Free energy landscapes of the May-Leonard model for various values of the diffusion constant D. The 
dynamics of the overall densities a, b, c is strongly confined to the invariant manifold of the well-mixed model, Eq. (17). To 
study the mechanisms underlying the distinct spatio-temporal patterns found in the spatial, stochast ic M ay-Leonard model 
we projected the probability for the overall densities onto the invariant manifold of the rate equations (17) for different values 



of the diffusion constant D. Color (gray scale) denotes the logarithmic probability to find the system globally in a specific 
state before reaching one of the absorbing states, such that red (medium gray) denotes a high probability, yellow (light gray) a 
medium probability and blue (dark gray) a low probability. The absorbing states themselves are not part of the statistics. For 
large D, no stable spatial structures can form and the dynamics corresponds to the well mixed case, Eqs. (17). As D becomes 



smaller than D « 9.5 • 10 - an attractor of the global dynamics emerges, effectively stabilizing the system against extinction. 
This attractor corresponds to planar, traveling waves with oscillating overall densities, and grows in radius with decreasing D 
due to a decreasing wavelength of the planar waves. For D < 3 • 10 -4 , a second attractor emerges, corresponding to rotating 
spirals. As a result of a decreasing wavelength of the spiral patterns the second attractor 's radius decreases with the diffusion 
constant, while the attractor corresponding to the traveling waves diminishes. Parameters were L — 60 and M — 8. Each plot 
was averaged over at least 100 realizations of the stochastic spatial dynamics. 



has been extensively studied in Refs. [34, 35 . Both pla- 
nar waves and rotating spirals are only met ast able, as 
stochastic fluctuations eventually lead to the annihilation 
of neighboring fronts and the dynamics will ultimately 
end in one of the three absorbing states which correspond 
to the extinction of two of the three species. The dynam- 
ics into the absorbing states has been found to be highly 
nontrivial, as the dynamical regimes described above lead 
to transitions into the absorbing states on different time 
scales. Furthermore, their statistical weight heavily de- 
pends on the diffusion coefficient D [30l [33 , 34, 38, Hi]. 



A. Global attractors and "free energy landscape" 
of the spatio-temporal dynamics 

To gain insight into the mechanisms responsible for 
these qualitatively different spatio-temporal patterns and 
how they determine the longevity of biodiversity in the 



population, we studied the global phase portrait of the 
dynamics. Figure [2] shows histograms for the overall con- 
centrations 

(a(t),b(t),c(t)) = J (a(x,£),6(x,£),c(x,£))d 2 x (18) 

of the three species on the invariant manifold of the rate 
equations, Eqs. (17a)-(17c). In detail, the negative log- 



arithm of the probability P(a, 6, c) to find the system in 
a specific global state (a, 5, c) before reaching one of the 
absorbing states is projected onto the invariant manifold: 



J 7 (a, 6, 



In P(a, 6, c), 



(19) 



The quantity T hence gives the logarithmic density of 
global trajectories in phase space and it can be consid- 
ered as an effective potential in the following sense: When 
instead of the mean-field reaction term, as given by the 
right hand side of Eq. (15), one uses J 7 as a "renormal- 



ized" potential for a Ginzburg-Landau theory one obtains 



7 



a reaction-diffusion equation which gives a good descrip- 
tion of the spatio-temporal dynamics. Long-lived spatio- 
temporal patterns correspond to regions of high proba- 
bility on the manifold, and are termed global attractors 
of the spatio-temporal dynamics in the following. Equiv- 
alently, adapting terminology from statistical mechanics, 
these attractors may be viewed as minima of the "free 
energy landscape" T . Of course, it is to be understood 
that these attractors are only metastable, i.e. while the 
system spends a long time in these states, ultimately de- 
mographic fluctuations will drive the system into one of 
the absorbing states. Intuitively, one may visualize those 
fluctuations as driving the escape of the dynamics from 
the minima in the "free energy landscape" into one of 
the absorbing states. We will later see that Qualitative 
changes in the shape of these minima correspond to tran- 
sitions in the nature of the dynamic processes leading 
to extinction and the mean times to extinction. These 
transitions should not be confused with non-equilibrium 
phase transitions. Rather they are to be considered as 
bifurcations in the nonlinear dynamics. 

Figure [2] shows that the shape of these global attrac- 
tors strongly depends on the magnitude of the diffusion 
constant D. For D > 10 -3 , there are no attractors other 
than the regions in the immediate vicinity of the three ab- 
sorbing states. All trajectories describing the global dy- 
namics quickly leave the unstable fixed point (a*,6*,c*) 
and approach the boundaries of the invariant manifold. 
Therefore, the probability is highest in the center (be- 
cause the dynamics starts there) and at the boundaries. 
In this regime, the system can be considered as well- 
mixed. The heteroclinic orbits in the global phase por- 
trait then correspond to spatially uniform oscillations be- 
tween states where one of the three species dominates; cf. 
Fig. [TJa). With decreasing diffusion constant D the na- 
ture of the global attractor changes qualitatively. Start- 
ing at the center of the manifold, the free energy develops 
a distinct local minimum which then evolves into a tri- 
angular shaped closed region; see Fig.[2{a)-(c). In other 
words, the phase portrait of the global population dy- 
namics changes from an unstable fixed point with hete- 
roclinic orbits to a pronounced limit cylce. Inspecting the 
spatio-temporal patterns as obtained from our stochastic 
simulations, we find that this limit cycle of the global 
dynamics corresponds to planar traveling waves; see also 
Fig.JTJb). The triangular shape is the result of oscillating 
overall densities. In the following we will refer to this par- 
ticular limit cycle as the wave attractor. Further lowering 
the diffusion constant, a second global attractor emerges 
as a smaller triangle inside the triangle corresponding 
to the planar waves, cf. Fig. [2jd)-(f). The inner trian- 
gular attractor corresponds to rotating spirals and will 
henceforth be denotes as the spiral attractor. We find 
that in this regime of diffusion constants the two attrac- 
tors coexist, meaning that we observe both planar waves 
and rotating spirals. Both processes may even be found 
within the same realization. With even further decreas- 
ing the diffusion constant, the weight of the triangular- 



shaped attractor corresponding to planar waves decreases 
and the attractor eventually disappears completely. As a 
consequence, the attractor of spiral waves gains weight, 
such that the dynamics at low mobilities is dominated by 
spiral waves. Taken together, we find that the phase por- 
trait of the global dynamics changes qualitatively upon 
decreasing the diffusion constant, and that those qual- 
itative changes have a one-to-one correspondence with 
distinct spatio-temporal patterns in the population dy- 
namics. As a consequence, the free energy landscape on 
the invariant manifold can be taken as a fingerprint of 
the spatio-temporal dynamics. We will use it in the fol- 
lowing to identify transitions between different patterns 
and analyze the ensuing changes in the dependence of 
the extinction times on system size. 

B. Pattern selection and extinction times 

The attractors in Fig. [2] show a triangular symmetry. 
A reduced representation for the global dynamics on the 
invariant manifold can therefore be obtained in terms of 
a properly defined radial variable. A convenient choice is 
the Lyapunov function 

a be . , 

C^——-, (20) 

(a + b + c) 

evaluated with the global concentrations a, 6, c. It mea- 
sures the distance of a global state to the boundaries 
of the invariant manifold and is approximately constant 
along the attractor for the planar waves. Figure [3] shows 
the effective free energy J^a, 6, c) as a function of the 
Lyapunov function and the diffusion constant. One easily 
identifies two threshold values of the diffusion constant 
where there are qualitative changes in the free energy 
landscape. We recover a threshold value D\ « 9 • 10 -4 
marking a transition from a well-mixed dynamics to a 
dynamics with spatio-temporal patterns [33j [34] . How- 
ever, the range of patterns is much richer than previously 
noted. Actually, the first threshold D\ marks a transition 
from spatially uniform oscillations between states domi- 
nated by a single species to planar waves where the three 
species cyclically chase each other. Note that the global 
oscillations still form part of the dynamics, albeit with a 
lower probability. Upon lowering the diffusion constant 
below a second threshold value, D 2 ~ 4.5 • 10 -4 , the his- 
togram of system trajectories becomes bimodal with a 
second metastable attractor emerging which is located 
close to the center of the invariant manifold. It can 
hence be identified with the inner, triangular attractor 
on the invariant manifold. As discussed before, this sec- 
ond attractor corresponds to rotating spiral waves. The 
coexistence of two attractors in this regime of mobili- 
ties means that both, planar waves and rotating spirals, 
are observed. Depending on the choice of initial condi- 
tions, the dynamics may at first end up in either one 
of the two attractors. Due to stochastic fluctuations it 
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FIG. 3. (Color online) Free energy landscape of the global 
dynamics of the two-dimensional May-Leonard system. The 
value of the Lyapunov function £ is a measure of how close 
a specific state is to the boundary of the invariant manifold 
(C = 0). Color (gray scale) denotes the logarithmic prob- 
ability to find the system at a specific value of £, i.e. the 
free energy formally defined in Eq. (19). Red (medium gray) 
signifies small values of the free energy (minima of the poten- 
tial) and thereby an attractor of the global dynamics. Yellow 
(light gray) denotes intermediate values, and blue large values 
of the free energy. The free energy landscape changes quali- 
tatively at two threshold values for the diffusion constant D. 
For large values of D the effective free energy has minima in 
the center (C = 0.037), where the dynamics starts, and at 
the boundaries of the invariant manifold (C = 0). At a first 
threshold D\ w 9 • 10 -4 an attractor emerges, which moves 
away from the reactive fixed point (C = 0.037) with decreas- 
ing values of D. Below a second threshold, D « 4.5 • 10 -4 , 
a second attractor emerges near the reactive fixed point, co- 
existing with the first one. For even smaller mobilities the 
dynamics is solely determined by the attractor near the reac- 
tive fixed point. Comparing with our simulations we find that 
these attractors correspond to global oscillations (heteroclinic 
orbits), planar waves and rotating spirals, respectively. The 
stochastic simulations were performed on a square lattice of 
linear size L = 60 and with a carrying capacity M — 8 for 
each site. For each values of D, the histogram was averaged 
over at least 100 realizations. 



may, however, from time to time switch between the two 
attractors akin to thermal fluctuations causing rare tran- 
sitions between different potential minima. With further 
decreasing D we observe that the metastable attractor 
corresponding to planar waves dissolves and only the at- 
tractor corresponding to rotating spirals remains. 

To further scrutinize the effect of these spatio-temporal 
patterns and the ensuing metastable global attractors on 
the system's dynamics, we analyzed the mean first pas- 
sage time into the absorbing states as a function of D, 




D 



FIG. 4. Mean lifetimes (dark dots) and coefficient of variation 
(gray squares) of the two-dimensional May-Leonard system as 
a function of the diffusion constant D. The mean lifetime in- 
creases abruptly at a first threshold value D± (indicated by a 
dashed line) where planar waves form. After passing through 
a maximum the lifetime decreases again. This is due to pla- 
nar waves which become unstable with decreasing correlation 
length. At the second threshold D2 (dashed line) rotating 
spirals become possible. With further decreasing diffusion 
constant the mean lifetime is asymptotically described by a 
power law D~ 2 . The coefficient of variation is a dimensionless 
measure for the dispersion of the probability distribution of T. 
The dispersion near the upper threshold D\ becomes large, 
i.e. we observe dynamical regimes on a variety of different 
time scales. 



see Fig. [4] We find that the mean time to extinction 
increases abruptly at D\, where the global attractor of 
planar waves emerges. After passing a peak value the 
mean lifetime then decreases again, as the wavelength 
of the planar waves becomes smaller. Then, as a result, 
these waves become more prone to fluctuations, and the 
rate of domain annihilation increases. Finally, below D2 
the lifetime increases again, which we attribute to the 
emergence of stable spiral waves. For small values of D 
the mean lifetime follows a power law (T) oc D~ 2 . This 
dependence can be understood by a simple scaling argu- 
ment: Since spirals annihilate pairwise as they meet, the 
mean lifetime should scale quadratically with their num- 
ber, (T) oc (n S pi ra i s ) 2 . The number of spirals in the sys- 
tem scales with their wavelength as ^spirals °c 

A" 2 . With 

A oc \/D we then infer that the mean lifetime scales as 
(T) oc D~ 2 , which is in good agreement with our numer- 
ical results. 

Figure [4] also shows the coefficient of variation, defined 
as the standard deviation divided by the mean, 



c v ee ^((T-(T)) 2 ) J (T) 



(21) 



It gives a dimensionless measure for the dispersion of the 
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probability distribution of T. We find that the disper- 
sion increases drastically right at the threshold D\. In 
this regime the standard deviation is much larger than 
the mean. From the spatio-temporal dynamics observed 
in our simulations we infer that this is due to the fact 
that there are several distinct dynamic processes driving 
the system towards an absorbing state and that these 
processes occur on greatly different time scales. There 
are rapid extinction processes, where - after a short tran- 
sient - domains in a planar wave are aligned in a non- 
cyclic order and thus immediately annihilate. We also 
find a process, where the global dynamics performs het- 
eroclinic orbits. Last, one observes metastable planar 
waves, cf. Fig. [I] Note that although the planar waves 
process is metastable it does not necessary mean that 
it dominates the long time properties of the system. In 
Ref. [38] it has been shown for the one dimensional model 
that probability of extinction scales differently with the 
system size for these two processes. In particular, one 
observes a crossover, such that for small systems planar 
waves determine the long time tails, while for very large 
systems global heteroclinic orbits are responsible for the 
longest living states. The relative weight of these pro- 
cesses depends on the diffusion coefficient D. As we have 
already learned from the above analyses, below the lower 
threshold value, D 2l there are also spiral waves emerg- 
ing. With decreasing D, spirals become the dominant 
patterns while all the other dynamic processes become 
less and less probable. As a result, the mean time to 
extinction is dominated by an escape out of the spiral 
attractor. The dispersion therefore decreases again. 

The probability distributions of first passage times of 
the above dynamical processes leading into the absorbing 
states show significantly different scaling behavior with 
the system size, cf. Ref. [38]. From an evolutionary per- 
spective, the tails of these distributions are most rele- 
vant because they correspond to rare, but extremely long- 
living communities maintaining biodiversity. The reason 
for their relevance is that the probability to observe a 
short-living (transient) ecosystem in nature is much lower 
than the probability to observe an ecosystem which per- 
sists for a long time. In Ref. [38] two of the authors 
showed that the tail of the distributions of first passage 
times of heteroclinic orbits scale like exp(T/7V), while 
for traveling waves the tail scales like exp [T/(ln7V) 3 ]. 
As a consequence, there is a crossover in the tail of the 
overall distribution of first passage times. Interestingly, 
while for small systems the long time dynamics is dom- 
inated by traveling waves, for large systems it is domi- 
nated by heteroclinic orbits. Although the computation 
of the distribution of first passage times is not feasible 
in two dimensions, we expect that similar arguments will 
hold here, as well. 

As shown in Ref. [33, 34 , there is a transition from a 
spatially uniform dynamics reminiscent of a well-mixed 
system to a dynamics dominated by spatio-temporal pat- 
terns when the wavelength of the pattern exceeds the 
system size. Following the classical theory of front prop- 



agation into unstable states [64], the wavelength of the 
traveling and spiral waves can be determined using the 
complex Ginzburg-Landau equation ( 15 ) [33j [36] 



A 



27TC3 



(1 - V 1 + c . 



(22) 



Due to the difference in geometry between planar and 
spiral waves this implies two distinct thresholds, D\ and 
D2. For planar waves on a square lattice we simply 
have the condition that the wavelength equals the sys- 
tem size, A(.Di) = 1. In Refs. [33] [34] it was found 
that the calculated wavelength deviates by a constant 
factor of 1.6 from the numerical value of the wavelength. 
This rescaling factor accounts for the renormalization of 
the reaction term due to spatio-temporal correlations, 
as captured by the global attractors. Using this rescal- 
ing factor we find a threshold value D\ ~ 7.6 • 10 -4 , 
in good agreement with the numerically found value, 
Di«9 - 10- 4 , cf. Fig. d The very same threshold is 
also found in the one dimensional model [38 . There, 
planar waves are the only possible spatial pattern and 
the threshold stems from their wavelength outgrowing 
the system size. Remarkably, the numerical values for 
D\ coincide in both, one and two spatial dimensions, as 
the complex Ginzburg-Landau equation predicts equal 
wavelengths for both cases, see further below. Since spi- 
rals always arise as pairs of anti-rotating spirals, stable 
pairs are possible, as long as the minimum distance d m i n 
between two vertices of the spirals is smaller than half 
of the system size. In other words, the threshold D 2 is 
given by d rn i n {D 2 ) = 1/2. To obey geometric constrains 
dictated by the periodic boundary conditions and the 
spirals' wavelength, the minimum distance of two anti- 
rotating spirals is d m i n = 2/3A(D), cf. also Fig. [ljd) . 
This implies a threshold value of A = 3/4 which is close 
to A ~ 0.8 obtained numerically in Ref. [34 . Hence, 
from 2/3A(L> 2 ) = 1/2 we obtain D 2 ~ 4.3 • 1(T 4 , in good 
agreement with the numerical results shown in Fig. [3] 

In the following we provide a scaling argument giving 
the scaling of the size of the wave attractor with the mo- 
bility D. In the intermediate regime between the two 
threshold values of D the wavelength of the planar wave 
patterns is of the same order as the system size. Hence, 
in this regime the finite spatial extension of the system 
is important. In our case, periodic boundary conditions 
allow stationary domain profiles only for certain values of 
the wave length, A = 1/n, n = 1, 2, . . . [65^ If the wave- 
length does not match any of these values, we observe 
oscillations in the overall concentrations, corresponding 
to the triangular attractor in Fig. [2] In the intermediate 
regime, D\ > D > D 2 , where A is slightly smaller than 
1, two domains take the characteristic domain size dic- 
tated by D and the third domain occupies the rest of the 
system. We now employ these intuitive observations to 
obtain the scaling of the wave attractor. Our numerical 
simulations reveal that the radius of the wave attractor 
increases with D according to a power law with an expo- 
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FIG. 5. (Color online) Illustration of the instability of wave 
fronts in the cyclic Lotka-Volterra system. The initial con- 
dition at t = was chosen as three domains of equal size in 
cyclic order. The pictures show snapshots at times t = 135, 
180, and 300. Color (gray scale) denotes species concentra- 
tions as described in Fig. [I] Parameters were D — 10 -4 , 
M = 8 and L = 80. 



nent of approximately 0.9, meaning that the correspond- 
ing values of the Lyapunov function increases with this 
exponent. If the system size is not a multiple of A ~ D 1 / 2 
one of the three domains will be of different wavelength. 
Assuming the concentration of empty sites is independent 
of D we set without loss of generality a = b ~ D~ x l 2 and 
c ~ 1 — 2D -1 / 2 . Upon inserting these concentrations into 
the Lyapunov function, we obtain C~ D- 2L> 3 / 2 , which 
is to leading order in agreement with our results [66] . 

Summarizing, we find that the spatio-temporal dynam- 
ics changes qualitatively at certain threshold values of the 
diffusion constant. These changes are finite size effects 
in the sense that they arise as a result of the comparison 
of certain length scales. We use the term "transition" 
for this behavior in the sense that macroscopic proper- 
ties of the system change qualitatively and abruptly at 
these threshold values. This is particularly evident in the 
mean first passage times to extinction. In the language 
of nonlinear dynamics the system undergoes bifurcations 
as a function of the mobility. 



III. THE CYCLIC LOTKA-VOLTERRA LIMIT 

In the limit a — » 0, ji —> only reactions remain, 
where the replication of predators does not require the 
availability of empty spaces. The resulting model is then 
of the Lotka-Volterra type [43] 04] , and characterized by 
a reduced set of chemical reactions: 



AB A AA, BC A BB, CA^CC. 



(23) 



This model is often referred to as the three-species Lotka- 
Volterra model. Although at a first glance there are 
no dramatic differences to the May-Leonard reactions, 
Eq. (16), the ensuing nonlinear dynamics is vastly differ- 



ent. The deterministic rate equations read 

d t a = va(b — c) , 
d t b = vb(c — a) , 
dtc = vc(a — b) . 



Without loss of generality, we also fix the normalization 
of total concentrations: a + b + c = 1. The nonlinear 
dynamics of the well-mixed cyclic Lotka-Volterra model 
again exhibits the same absorbing fixed points as the gen- 
eral model ([3]). The reactive fixed point is now given by 



(25) 



It is, however, neutrally stable as the real parts of the 
eigenvalues, Eq. ([H]), are zero. In fact, C = for any a, 
b and c, such that starting from any point on the phase 
plane, the trajectories form neutrally stable cycles. 

Similar to the May-Leonard model, species' mobility 
drastically alters the system's collective dynamics [26] . 
However, the ensuing spatio-temporal dynamics of the 
cyclic Lotka-Volterra and the May-Leonard model differ 
qualitatively. This behavior can be understood upon con- 
sidering the dynamics of domain boundaries separating 
different species. In the May-Leonard model the separa- 
tion of selection and reproduction processes is counter- 
acting the roughening of these domain boundaries due 
to stochastic fluctuations: If a species from one domain 
enters the other species' domain, it first creates empty 
sites. Since these empty sites are occupied with a higher 
probability by offsprings of individuals from the invaded 
species rather than by invaders, the invasion process is 
unlikely to be successful. This stabilizes spatially sepa- 
rated domains in the May-Leonard model. In contrast, 
in the cyclic Lotka-Volterra model an invader directly 
replaces the invaded species such that it has a higher 
probability of success. As illustrated in Fig. [5] this leads 
to a roughening instability of planar wave fronts. How- 
ever, this does not imply the total loss of any spatial 
correlations. To the contrary, there are still strong corre- 
lations and they play a fundamental, yet subtle, role in 
the spatio-temporal dynamics and the processes leading 
to the extinction of all but one species. 

In our simulations we observe different dynamic pro- 
cesses depending on the mobility. For large diffusion con- 
stants, where the system can be considered well- mixed, 
we recover the homogeneous oscillations as predicted by 
the rate equations ( 24 ) . We still find homogeneous oscil- 



(24) 



lations if the mobility is decreased. However, as we will 
see below, these systemwide oscillations are of entirely 
different nature as the neutrally stable orbits found in 
the well- mixed system. For even lower mobility we finally 
find a seemingly random appearance and dissolution of 
spatial clusters. These clusters are convectively unstable 
spiral waves, which, due to a roughening transition as- 
sociated with an Eckhaus instability, appear, move and 
annihilate. 



A. Extinction times and extinction probabilities 

As discussed previously [26] 33, 34!, a convenient mea- 
sure to characterize the stability of the system is the 



11 



probability P ext that the system has reached an absorb- 
ing state within a time proportional to the system size 
N. The simulations for our model reproduce the results 
found in Ref. [26 j. For large D our result coincides with 
the analytically obtained value found for the non-spatial 
system [17] . see Fig. 6^a). For low mobilities the extinc- 
tion probability is close to zero. Hence, the system is 
in a metastable state with extinction times scaling expo- 
nentially in the system size, cf. Fig. [6jb). At a thresh- 
old value D c « 3 • 10 -3 , there is a sharp transition to 
P ext ~ 1 indicating that extinction times scale logarith- 
mically in the system size N. Indeed, Fig. |6jc) indicates 
that the scaling of mean first passage times is sublin- 
ear. For even larger values of the diffusion constant, the 
extinction probability decreases again until it reaches a 
value of 0.8 in the well-mixed limit. Here, extinction 
times scale linearly, as demonstrated by Fig.[6jc). Hence, 
the global dynamics is characterized by the escape out of 
a neutrally stable state. We conclude that spatial cor- 
relations increase the system's stability for small D and 
destabilize it above a threshold value D r . 



B. Global attractors and free energy landscapes 

As for the May-Leonard model we now employ a study 
of the global phase portraits to gain a deeper understand- 
ing of the ambiguous impact of spatial structures on the 
longevity of biodiversity. The existence of metastable 
states below a certain mobility threshold, suggested by 
the scaling of extinction times with system size, is sup- 
ported by histograms of the global dynamics. Figure [7] 
shows the free energy landscape projected onto the in- 
variant manifold of Eqs. (24). For very small values of D 



we find an attracting region in the center of the simplex. 
This attractor corresponds to convectively unstable spi- 
rals. As mentioned before, smooth domain borders are 
subject to roughening and therefore become unstable in 
the Lotka-Volterra model. However, while spatial pat- 
terns can not be maintained, strong correlations exist 
and effectively render the global dynamics metastable. 
With increasing values of D we observe that the trajec- 
tories describing the overall dynamics of the system are 
attracted towards a limit cycle, which grows in radius and 
eventually reaches the boundary of the invariant mani- 
fold. This limit cycle corresponds to systemwide oscil- 
lations. As these oscillations are linked to a metastable 
attractor, they are much more long-lived compared to 
the neutrally stable oscillations found in the well-mixed 
case, i.e. their mean lifetime scales exponentially with 
the system size. At some threshold value of the diffusion 
constant the attractor coincides with the boundary of the 
simplex. Then, the absorbing states are embedded within 
the limit cycle. As a consequence, the global dynamics is 
effectively attracted towards the boundaries of the phase 
plane once it reaches the limit cycle's basin of attraction, 
and therefore rapidly reaches one of the absorbing fixed 
points. The global dynamics is therefore effectively het- 



eroclinic and approaches the absorbing states exponen- 
tially fast. Hence, in this regime spatial structure desta- 
bilize the system, which explains the sub linear scaling of 
extinction times as shown in Fig.[6ja) and (c). For large 
D, the attractor lies outside of the simplex, such that the 
global dynamics on the simplex is essentially governed by 
neutrally stable orbits. 

Figure [8] illustrates the different dynamical regimes by 
means of the effective free energy T for the cyclic Lotka- 
Volterra model. We identify three distinct regime, which 
are characterized by the shape of the effective free energy. 
For the well-mixed system, D > D c « 3- 10 -3 , the poten- 
tial is flat and the global dynamics is neutrally stable as 
predicted by the rate equations ( 24 ) . At D c an attractor 



emerges, which at this point coincides with the bound- 
aries of the simplex (C = 0). With decreasing values of D 
the attractor is located at increasingly large values of the 
Lyapunov function until at D « 4 • 10 -4 it coincides with 



the reactive fixed point of the rate equations (24). Com- 



paring with our simulations we therefore identify three 
regimes: neutrally stable orbits, metastable systemwide 
oscillations, and convectively unstable spirals. In conclu- 
sion, the behavior of the attractors of the global dynamics 
provides an intuitive explanation for the observed tran- 
sitions in the extinction probabilities. 



IV. THE INTERMEDIATE REGIME 

While the previous sections considered important lim- 
iting cases of the reactions (J3|, we now study the gen- 
eral case with cr, /i, v ^ 0. We will use the following 
parametrization, which allows to tune the relative weight 
of Lotka-Volterra type reactions and May-Leonard type 
reactions: 



i/(s) = s, 
ee 1, 
a(s) ee 1 - 



(26) 



Here the parameter s is the fraction of Lotka-Volterra 
type reactions, and is varied between and 1. This choice 
of parametrization has two important properties: First, 
it conserves the limits discussed in the previous sections 
and makes them comparable. In the Lotka-Volterra limit 
and in the May-Leonard limit per time step each individ- 
ual performs, on average, one active selection process or 
passive process, respectively. This holds for any value 
of s. Second, our simulations show that the correlation 
length of species concentrations stays approximately con- 
stant when changing s (data not shown here). This is 
because in our parametrization we fix the relevant time 
scale, and thereby by dimensional analysis, for a given 
mobility, the correlation length. 

Figure [9] shows that with increasing values of s spiral 
patterns become convectively unstable, i.e. the vertices 
start to move and annihilate. The destabilizing effect 
of Lotka-Volterra reactions on spiral patterns can also be 
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FIG. 6. (a) Probability that the system with Lotka-Volterra reactions reaches an absorbing state before t — N. We observe a 
sharp transition from survival to extinction at D « 3 • 10 -3 . In the well mixed limit the extinction probability converges to a 
finite value of 0.8 (M = 8, L = 60). (b) and (c) show the scaling of the mean first passage time into any of the absorbing states 
with the system size N. Two phases can be identified. For D < 3 • 10 -3 the scaling becomes exponential, hinting at an escape 
process from a metastable state. In the well mixed case (D — 100) the scaling is linear, in agreement with the escape out of a 
neutrally stable state. In the intermediary regime our results are in agreement with both a logarithmic and a linear scaling. 
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FIG. 7. (Color online) Probability to find the system glob- 
ally in a specific state before reaching one of the absorbing 
states. Red (medium gray) denotes a high probability, yel- 
low (light gray) a medium probability, and blue (dark gray) a 
low probability. One observes the emergence of an attracting 
limit cycle of the global dynamics. The attract or grows in 
radius with increasing D and eventually reaches the bound- 
aries of the invariant manifold. For even larger values of the 
diffusion constant the attractor lies outside of the invariant 
manifold and the global dynamics is essentially neutral. The 
histograms were sampled over 100 trajectories until T — N . 
Parameters were M — 8 and L — 60. 



species is present at approximately equal concentrations, 
and therefore indicate the position of spiral vertices. In 
Fig. |9|d)-(f) black dots correspond to positions, where 
this absolute value is smaller than 0.13 [67]. We thus 
infer that the spirals become unstable with increasing 
s. Indeed, the complex Ginzburg-Landau equation (15) 



predicts an Eckhaus instability, implying that the spi- 
rals vertices become convectively unstable [26 , i.e. they 
move, annihilate and appear again constantly above a 
certain value of s. To determine this value we follow the 
steps given in Ref. [55 , where the stability of planar wave 
solutions was studied. The waves are stable, as long as 
the generalized Eckhaus criterion, 



1 2 (1 + C ^% 
1 1-Q 2 >U ' 

holds, where Q is the selected wave vector, 



D 

ci(s) 



yi+c 3 (s) 2 



(27) 



(28) 



Inserting cs(s) and solving for s we find a critical value 
of se ~ 0.32. The breakdown of stable spatial struc- 
tures as the result of a roughening transition is indeed 
confirmed by our numerical simulations. In contrast to 
the transitions in D, the Eckhaus instability is indepen- 
dent of the size of spatial patterns and it can therefore 
be considered a transition in the strict thermodynamic 
sense. It has significant, yet ambiguous, implications for 
the stability of biodiversity, as will be discussed in the 
following. 



visualized by considering the absolute value of the coordi- 
nates defined in Eq. ([7]), |y(a, 6, c)\. It gives a measure of 
how far the system is locally away from the reactive fixed 
point. Low values of |y| correspond to a locus where each 



A. Extinction times 



Figure 10 shows the mean first passage time to one of 
the absorbing states as a function of D and s. The color 



13 



123456789 10 11 




10" 1 10" 2 10" 3 10" 4 



FIG. 8. (Color online) The global phase portrait of the Lotka- 
Volterra system. For each diffusion constant D we plot (in col- 
orcode or gray scale) the probability to find the system before 
t = N at a specific value of the Lyapunov function C = abc. 
Red (medium gray) denotes a high probability, yellow (light 
gray) a medium probability, and blue (dark gray) a low prob- 
ability. Three dynamical regimes corresponding to neutrally 
stable orbits, systemwide oscillations and convectively unsta- 
ble spirals can be identified and linked to their corresponding 
attractors of the global dynamics. Note that the attractor 
vanishes at D « 3 • 10 -3 corresponding to the abrupt increase 
of extinction probabilities in Figure [6] The simulations were 
done with M = 8 and L = 60, and, for numerical reasons, 
stopped at T = N. For each of roughly 20 data points in D 
we averaged over approximately 100 trajectories. 



code as indicated in the figure is chosen such that red 
corresponds to large and blue to short extinction times. 
Dark red indicates the longest time simulated, t — 10 T . 
The limits s = and s = 1 correspond to the May- 
Leonard and Lotka-Volterra models, respectively. Vary- 
ing 8, however, does not simply interpolate between these 
two limits, but leads to a rather rich and complex dynam- 
ics. In particular, there is a local maximum in the mean 
extinction time for finite values of s below se- We in- 
fer from our simulations that this maximum is linked to 
the emergence of planar traveling waves. In contrast to 
the May-Leonard model (s = 0), planar waves seem to 
be increasingly important in this regime: They dominate 
the dynamics for a rather broad range in the diffusion 
coefficient. Moreover, they seem to be much more stable 
as compared to the May-Leonard case, which can be seen 
by comparing Fig. [4] and Fig. [TUJ While the exact reason 
for this remains unclear, the stabilization of planar waves 
seems to be related to a change in the wavelength and 
thereby a reduction in the oscillations of the global con- 
centrations. This can be inferred from the global phase 
portraits, as discussed below. For small D, we again find 



8=0 8=0.5 s=l 




FIG. 9. (Color online) (a)-(c) Snapshots of the spatial distri- 
bution of species for different values of the fraction of Lotka- 
Volterra reactions s indicated in the graph. Color denotes 
(gray scale) the concentration of the species A, B and C, as 
described in Fig. ^ With increasing s spirals become convec- 
tively unstable, i.e. they move, annihilate and then appear 
again, (d)-(f) To illustrate the destabilization of spiral waves 
we computed for each lattice site the distance from the re- 
active fixed point |y(a, 6, c)|. Dark points show sites where 
|y(a, b, c)| is below a certain threshold, thereby indicating the 
position of spiral vertices. Parameters were D — 10 -4 (corre- 
sponding to the regime, where spirals and waves are possible 
in the May-Leonard model), M — 8 and L — 80. 



metastable rotating spirals. For the well-mixed system 
we find short first passage times. The concentrations 
there perform homogeneous oscillation, which we identi- 
fied with heteroclinic orbits of the global trajectories for 
the May-Leonard case 8 = 0. These orbits cover a broad 
parameter regime. In particular, they also arise for values 
of 8, where most of the reactions are of Lotka-Volterra 
type. The reason for this can be inferred from the stabil- 
ity of the reactive fixed point of the rate equaltions Q. 
The corresponding eigenvalues ([6| retain a non vanishing 
positive real part. The trajectories of the global dynam- 
ics are therefore driven to the vicinities of the absorbing 
points exponentially fast. As a result, even for s « 0.9 
the global dynamics is determined by a tiny fraction of 
May-Leonard reactions. 

The roughening transition is complicated by threshold 
values in D, corresponding to the onset of planar waves 
and spirals, and the dissolution of the former. These 
threshold values take the same values as in the limiting 
case of only May-Leonard reactions. As the value of s 
exceeds the roughening transition (Eckhaus instability) 
we observe a sharp transition between long extinction 
times for small values of D and short extinction times for 
large values of D. In the latter regime, spirals and pla- 
nar waves become convectively unstable as predicted by 
the complex Ginzburg-Landau equation ( fl5| ) . For spiral 
waves this is illustrated by Fig. [9] Nevertheless, strong 
correlations exist and mean times to extinction are large 
in this regime. From our simulations we infer that the 
dominant dynamic process in this regime can be identi- 
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FIG. 10. (Color online) Average first passage time into any of 
the absorbing states as a function of the diffusion constant D 
and the fraction of Lotka-Volterra reactions s. Red (medium 
gray) denotes a large lifetime, yellow (light gray) a medium 
lifetime, and blue (dark gray) a small lifetime. For s = we 
obtain the mean lifetimes shown in Fig. [4] The dynamics is es- 
sentially governed by heteroclinic orbits, traveling waves and 
rotating spirals. With increasing s the planar waves become 
more and more stable and dominate the dynamics for a full 
order of magnitude in D. The prominence of traveling waves 
leads to a local maximum in the mean lifetimes. For even 
higher s the system undergoes an Eckhaus instability (the 
analytical value is denoted by a dashed line), where planar 
waves become unstable. The dynamics is roughly comparable 
to the heteroclinic orbits in the May-Leonard model. Neutral 
orbits are driven to the boundary of the invariant manifold 
by a limited fraction of May-Leonard reactions. For s = 1 we 
again recover the dynamics of the Lotka-Volterra model stud- 
ied in Section |III| For each of the approximately 400 data 
points averages were taken over about 100 trajectories. Due 
to numerical constrains simulations were stopped at T = 10 7 . 
Parameters were M = 8 and L = 60. 



fled as the convectively unstable spirals also found in the 
Lotka-Volterra limit. Note, however, that due to a trun- 
cation of simulation times not all details may be resolved 
in this regime. 



B. Effective free energy and Lotka-Volterra limit 

To study how the Lotka-Volterra limit is reached, we 
computed the effective free energy J 7 as a function of C. 
We focus on the case D = 3 • 10 -4 , which entails the 
regime of stable planar waves, cf. Fig. 10 In the May- 



Leonard model this corresponds to the regime shortly 
below the lower critical point in the diffusion constant, 
where the wave attractor and the spiral attractor coex- 



ist. Figure [TT] demonstrates that the observed changes 
in extinction times are related to the emergence, disap- 
pearance and changes in the characteristics of attractors 
of the global dynamics. The limit of the May-Leonard 
model (s = 0) was already discussed in Section [TTJ At- 
tractors for rotating spirals and planar waves are visible. 
When the fraction of Lotka-Volterra reactions is slightly 
increased the spiral attractor disappears while the wave 
attractor remains, see Fig. [Tljb) . The latter shrinks in 
size, hinting at an increasing wave length [Figs. [TT^c) -(e)]. 
The attractor then contracts towards the reactive fixed 
point [Figs. [OJd)-(f)]. In Fig. [To] this regime corresponds 
to the local maximum in extinction times. At the point, 
where the system undergoes an Eckhaus instability, the 
attractor dissolves [Fig. [TT^g) and (h)]. The dynamics is 
then dominated by global oscillations which are driven 
outward by a limited number of May-Leonard reactions 
[Fig.[]]Jh) and (i)]. This regime is therefore closely re- 
lated to the heteroclinic orbits found in the May-Leonard 
model. For an even larger fraction of Lotka-Volterra re- 
actions a new attracting limit cycle emerges, correspond- 
ing to the systemwide oscillations found in Section III 

[Fig.intD-a)]. 



The results are summarized in a free energy landscape 
as a function of s, cf. Fig. [12] For s = we find the 
attractors of the planar and spiral waves of the May- 
Leonard model, cf. Fig. [3] When the fraction s of Lotka- 
Volterra reactions is increased the attractor of planar 
waves shrinks to the center of the manifold. As a result, 
there are no oscillations in the overall densities, which is 
in contrast tothe May-Leonard model, where these oscil- 
lations stem from waves having a wavelength close but 
unequal to the system size. As a result of the lack of 
oscillations, planar waves become increasingly stable in 
this regime. At the Eckhaus instability, se, spatial pat- 
terns become unstable. The dynamics can then be best 
described as heteroclinic orbits. The system globally per- 
forms orbits, which are driven to the boundary of the 
manifold by the reactions of May-Leonard type. Hence, 
even a tiny fraction of May-Leonard reactions determines 
the global dynamics in this regime. This is not surprising, 
as the conservation law associated with the cyclic Lotka- 
Volterra model holds precisely only in the case s = 1. 
For values of s close to 1 reactions involving empty sites 
become unimportant. We then find the attractor corre- 
sponding to systemwide oscillations, cf. Fig. [8] Summa- 
rizing, in the model with direct and indirect competition 
we find a surprisingly rich variety of dynamic processes 
affecting the longevity of biodiversity in a much more 
complex manner than one would naively expect from an 
Eckhaus instability. In particular, we observed a local 
maximum in mean lifetimes if direct competition is weak 
but non vanishing. 
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FIG. 11. (Color online) Probability to find the system in a specific state before reaching one of the absorbing fixed points for 
fixed D = 3 • 10" 4 . The histo gram is projected onto the invariant manifold of the rate equations Q. We varied the fraction 
of Lotka-Voltera reactions from s — (top left) to s = 1 (bottom right). Starting at the classic May-Leonard model (s = 0), 
where attractors for planar waves and rotating spirals can be identified, the attractor for the spirals disappears with growing 
s. The remaining attractor contracts to the reactive fixed point and also, when the system undergoes an Eckhaus instability, 
dissolves. For an even larger fraction of Lotka-Volterra reactions the global dynamics is driven outward by a limited fraction 
of May-Leonard reactions and is comparable to the heteroclinic orbits found in the May-Leonard model. When a majority of 
the reactions are of Lotka-Volterra type, i.e. for s not much smaller than 1, we again observe the emergence of an attracting 
limit cycle corresponding to systemwide oscillations. Parameters where M — 8 and L — 60. 



16 



10 



12 



T{C) 



systemwide 
oscillations 




0.1 0.2 03 0.4 0.5 0.6 0.7 OA 

V7 



0.9 1 



FIG. 12. (Color online) To study the most intriguing line of 
Fig. 
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D = 3- 10 -4 , in more detail we computed the effective 
free energy T at specific values of the Lyapunov function C for 
different values of s. Red (dark gray) denotes minima of the 
potential and thereby attractors of the global dynamics. Yel- 
low (light gray) signifies intermediate values and blue (dark 
gray) large values of the effective free energy. We identify 
several regimes depending on the relative strength s of the 
different types of competition. For s = 0we recover the coex- 
isting wave and spiral attractors of the May-Leonard model. 
With increasing values of s only the wave attractor remains 
and approaches the reactive fixed point of the global dynam- 
ics (C = 0.037). At the Eckhaus instability (dashed line) the 
wave attractor dissolves. Instead, an attractor correspond- 
ing to global heteroclinic orbits emerges. Only, when almost 
all reactions are of Lotka-Volterra type an attractor close to 
the reactive fixed point emerges. The latter corresponds to 
the limit cycle found in the cyclic Lotka-Volterra model, see 
Sec. |III| Comparing with our simulations we find that these 
attractors are linked to rotating spirals, planar waves, global, 
heteroclinic orbits and systemwide oscillations. Simulation 
parameters were M = 8 and L = 60. 



V. CONCLUSION 

We studied the population dynamics of three-species 
models where species interact with each other cycli- 
cally through both direct and indirect competition. In 
the limiting cases of only direct or indirect competi- 
tion our model reduces to the cyclic Lotka-Volterra or 
May-Leonard model, respectively. For a well-mixed sys- 
tem, the nonlinear dynamics of these models differs sig- 
nificantly. While in both cases the trajectories lie on 
two-dimensional invariant manifolds comprising the ab- 
sorbing states of extinction, their phase portraits differ 
qualitatively. The dynamics of the cyclic Lotka-Volterra 
model is characterized by neutrally stable orbits. In con- 
trast, the dynamics of the May-Leonard model features 



an unstable fixed point in the center of the invariant man- 
ifold and heteroclinic orbits which approach the bound- 
aries of the invariant manifold and hence the absorbing 
states exponentially fast. In the spatial versions of these 
models, these attractors of the well-mixed system still 
act locally on each lattice site. However, if one views 
the spatially extended system as a set of interconnected 
local patches, the coupling between these patches due to 
diffusion may lead to qualitative changes in the type and 
stability of these attractors. 

Indeed, numerical simulations show that in spatially 
extended systems there is a rich diversity of spatio- 
temporal patterns depending on the systems' parame- 
ters. The goal of this work was to identify and char- 
acterize the dynamic processes responsible for the tran- 
sient maintenance of biodiversity and ultimately leading 
to extinction in the spatial models. To this end, we inves- 
tigated the phase portrait of the overall concentrations 
for the species comprising the system and analyzed the 
ensuing global attractors of the dynamics and how they 
are connected with the different types of spatio-temporal 
dynamics. Moreover, based on a statistical analysis of 
the system trajectories on the global phase portrait, we 
defined an effective free energy which gave us valuable in- 
formation about the scaling of extinction times with the 
system size. In particular, the minima of the free energy 
correspond to metastable dynamical processes. 

In the limit corresponding to the spatial May-Leonard 
model, the minima in the effective free energy landscape 
of the global phase portrait are linked to three distinct 
spatio-temporal patterns: (i) spatially homogeneous os- 
cillatory behavior, (ii) planar traveling waves and (iii) 
rotating spirals. Importantly, the characteristics of the 
global attractors change qualitatively at certain thresh- 
old values of the mobility. This means that the length 
scales associated with the spatial patterns changes, which 
affects their stability and thereby the probability to find 
the system in such a state. In particular, below an up- 
per threshold value of the mobility a triangular attractor 
corresponding to traveling waves emerges. This attrac- 
tor can be regarded as a limit cycle of the global dynam- 
ics. It grows in size for decreasing mobility, reflecting 
a decreasing wavelength. At a lower threshold value of 
the mobility a second limit cycle of the global concen- 
trations is found, which is located inside the attractor of 
the traveling waves. There, rotating spirals emerge. In 
this regime we observe the coexistence of two dynamic 
processes, planar waves and rotating spirals. Which of 
the two dynamic regimes is actually realized is deter- 
mined by stochasticity and the initial conditions. For 
even lower mobility, the attractor of the traveling waves 
dissolves, as the correlation length is too small compared 
to the system size to maintain planar domain borders. In 
this regime only the attractor of rotating spirals remains, 
which dominates the dynamics predominantly. 

As opposed to this behavior, in the limit of reactions 
of Lotka-Volterra type only, the system does not exhibit 
stable spatial patterns. However, there are strong spatial 
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correlations, and they manifest themselves in an attrac- 
tor of the global dynamics. This attractor has the form 
of a rounded triangle around the reactive fixed point and 
corresponds to a limit cycle of the global concentrations. 
The radius of the limit cycle grows with increasing mobil- 
ity: when the attractor, with increasing mobility, reaches 
the boundaries of the invariant manifold the dynamics 
passes from metastability (exponential scaling of extinc- 
tion times with the system size) to rapid extinction (sub- 
linear scaling of extinction times with the system size). 
For even larger mobilities, the radius of the global limit 
cycle outgrows the boundaries of the simplex. The mean 
time to extinction then scales linearly with the system 
size. Hence, in contrast to the May-Leonard model a sin- 
gle attractor here is responsible for three distinct scaling 
regimes. 

Finally, we found a remarkably complex behavior when 
varying the relative strengths of direct (Lotka-Volterra) 
and indirect (May-Leonard) competition. If direct com- 
petition is weak compared to indirect competition, pla- 
nar traveling waves are an increasingly important con- 
stituent of the extinction dynamics as compared to the 
May-Leonard case. These planar waves are very stable, 
leading to a local maximum of extinction times in the 
phase diagram. Simultaneously, we observe that in con- 
trast to the May-Leonard model rotating spirals do not 
form for intermediary mobilities. This is reflected in the 
dissolution of the corresponding attractor. At a specific 
fraction of Lotka-Volterra reactions the system undergoes 
an Eckhaus instability: traveling waves and rotating spi- 
rals become unstable. The Eckhaus instability manifests 
itself in the vanishing of the attractors of planar waves 
and rotating spirals. Beyond the Eckhaus instability, a 
new attractor emerges, corresponding to global oscilla- 
tory behavior for high mobility and convectively unsta- 
ble spirals for low mobilities. Summarizing, we find that 



the spatia-temporal dynamics of cyclic populations mod- 
els with both, direct and indirect competition is surpris- 
ingly rich and differs qualitatively from the cyclic Lotka- 
Volterra and May-Leonard models. We identified several 
threshold values of the mobility and the relative strength 
of the two types of competition. 

In conclusion, scaling of extinction times with the sys- 
tem size change abruptly at certain threshold values of 
the mobility and the relative strength of the two types 
of competition. We showed that the dynamic processes 
leading to the transient maintenance of biodiversity are 
linked to attractors of an effective free energy of the over- 
all concentrations. The characteristics of these attractors 
change upon certain threshold values, thereby giving in- 
sight into the mechanisms underlying these phase tran- 
sitions. By means of extensive numerical simulations we 
provide the complete phase diagrams, which are ratio- 
nalized by scaling arguments based on properties of the 
complex Ginzburg-Landau equation. 

We believe that the method of global phase portraits 
and the ensuing effective free energy landscapes (renor- 
malized reaction terms) might also give a deeper in- 
sight into the dynamics of spatial ecological models and 
reaction-diffusion systems in other fields of biology. In 
particular, further studies may apply this method to un- 
derstand epidemic models, asymmetric four species mod- 
els or more complex food webs. 
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